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C\l ' ABSTRACT 

In the Nice model, the late heavy bombardment (LHB) is related to an orbital instability of giant planets which causes a fast dynamical 

1 ■ dispersion of a transneptunian cometary disk. We study effects produced by these hypothetical cometary projectiles on main-belt 

asteroids. In particular, we want to check whether the observed collisional families provide a lower or an upper limit for the cometary 
flux during the LHB. 

We present an updated list of observed asteroid families as identified in the space of synthetic proper elements by the hierarchical 
clustering method, colour data, albedo data and dynamical considerations and we estimate their physical parameters. We selected 
12 families which may be related to the LHB according to their dynamical ages. We then used collisional models and N-body orbital 
simulations to gain insight into the long-term dynamical evolution of synthetic LHB families over 4 Gyr. We account for the mutual 
collisions between comets, main-belt asteroids, and family members, the physical disruptions of comets, the Yarkovsky/YORP drift 
in semimajor axis, chaotic diffusion in eccentricity/inclination, or possible perturbations by the giant-planet migration. 
Assuming a "standard" size-frequency distribution of primordial comets, we predict the number of families with parent -body sizes 
Z3pb > 200 km - created during the LHB and subsequent ^ 4 Gyr of collisional evolution - which seems consistent with observations. 
However, more than 100 asteroid families with D PB > 100 km should be created at the same time which are not observed. This 
discrepancy can be nevertheless explained by the following processes: i) asteroid families are efficiently destroyed by comminution 
(via collisional cascade), ii) disruptions of comets below some critical perihelion distance (q < 1.5 AU) are common. 
Given the freedom in the cometary-disruption law, we cannot provide stringent limits on the cometary flux, but we can conclude that 
the observed distribution of asteroid families does not contradict with a cometary LHB. 

i - ^ , Key words, celestial mechanics - minor planets, asteroids: general - comets: general - methods: numerical 

> ■ 

■ 1. Introduction asteroid samples heated by impact events, together with dynam- 
| ical modelling work, to suggest that the basin-forming portion 

The late heavy bombardment (LHB) is an important period in the of the LHB lasted from approximately 4.1-4.2 to 3.7-3.8 billion 

history of the solar system. It is often defined as the process that years ago on the Moon (Bogard 1995, 201 1, Swindle et al. 2009, 

\ made the huge but relatively young impact basins (a 300 km or Bottke et al. 2012, Norman & Nemchin 2012). 

■ larger diameter crater) on the Moon like Imbrium and Orientale. Tne so-called 'Nice model' provides a coherent explanation 
The sources and extent of the LHB, however, has been under- f me origin of the LHB as an impact spike or rather a "saw- 
going recent revisions. In the past, there were two end-member tootn » (Morbidelli et al. 2012). According to this model, the 
schools of thought describing the LHB. The first school argued bombardment was triggered by a late dynamical orbital insta- 

; that nearly all lunar basins, including the young ones, were made bility of the giant planets, in turn driven by the gravitational 

03 . by impacting planetesimals left over from terrestrial planet for- interactions between said planets and a massive transneptunian 

mation (Neukum et al. 2001, Hartmann et al. 2000, 2007; see disk of planetesimals (see Morbidelli 2010 for a review). In this 

Chapman et al. 2007 for a review). The second school argued scenario, three projectile populations contributed to the LHB: 

that most lunar basins were made during a spike of impacts that the comets from the original transneptunian disk (Gomes et al. 

took place near 3.9 Ga (e.g., Tera et al. 1974, Ryder et al. 2000). 2005), the asteroids from the main belt (Morbidelli et al. 2010) 

Recent studies, however, suggest that a compromise scenario and those from a putative extension of the main belt towards 

may be the best solution: the oldest basins were mainly made by Mars < inwards of its current inner edge (Bottke et al. 2012). The 

leftover planetesimals, while the last 12-15 or so lunar basins last could have been enough of a source for the LHB, as recorded 

were created by asteroids driven out of the primordial main belt in the lunar crater record (Bottke et al. 2012), while the asteroids 

by the effects of late giant-planet migration (Tsiganis et al. 2005, from the current main belt boundaries would have only been a 

Gomes et al. 2005, Minton & Malhotra 2009, Morbidelli et al. minor contributor (Morbidelli et al. 2010). 
2010, Marchi et al. 2012, Bottke et al. 2012). This would mean The Nice model, however, predicts a very intense cometary 

the LHB is limited in extent and does not encompass all lunar bombardment of which there seems to be no obvious traces on 

basins. If this view is correct, we can use studies of lunar and the Moon. In fact, given the expected total mass in the original 
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transneptunian disk (Gomes et al. 2005) and the size distribution 
of objects in this disk (Morbidelli et al. 2009), the Nice model 
predicts that about 5 x 10 4 km-size comets should have hit the 
Moon during the LHB. This would have formed 20 km craters 
with a surface density of 1.7 x 10~ 3 craters per km 2 . But the 
highest crater densities of 20 km craters on the lunar highlands 
is less than 2 x 10~ 4 (Strom et al. 2005). This discrepancy might 
be explained by a gross overestimate of the number of small bod- 
ies in the original transneptunian disk in Morbidelli et al. (2009). 
However, all impact clast analyses of samples associated to ma- 
jor LHB basins (Rring and Cohen 2002, Tagle 2005) show that 
also the major projectiles were not carbonaceous chondrites or 
similar primitive, comet-like objects. 

The lack of evidence of a cometary bombardment of the 
Moon can be considered as a fatal flaw in the Nice model. 
Curiously, however, in the outer solar system we see evidence 
of the cometary flux predicted by the Nice model. Such a flux is 
consistent with the number of impact basins on Iapetus (Charnoz 
et al. 2009), with the number and the size distribution of the 
irregular satellites of the giant planets (Nesvorny et al. 2007, 
Bottke et al. 2010) and of the Trojans of Jupiter (Morbidelli et 
al. 2005), as well as with the capture of D-type asteroids in the 
outer asteroid belt (Levison et al., 2009). Moreover, the Nice 
model cometary flux is required to explain the origin of the col- 
lisional break-up of the asteroid (153) Hilda in the 3/2 resonance 
with Jupiter (located at 4 AU, i.e. beyond the nominal outer 
border of the asteroid belt at 3.2 AU; Broz et al. 201 1). 

Missing signs of an intense cometary bombardment on the 
Moon and the evidence for a large cometary flux in the outer 
solar system suggest that the Nice model may be correct in its 
basic features, but most comets disintegrated as they penetrated 
deep into the inner solar system. 

To support or reject this possibility, this paper focusses on the 
main asteroid belt, looking for constraints on the flux of comets 
through this region at the time of the LHB. In particular we focus 
on old asteroid families, produced by the collisional break-up 
of large asteroids, which may date back at the LHB time. We 
provide a census of these families in Section [2] 

In Section [3] we construct a collisional model of the main 
belt population. We show that, on average, this population alone 
could not have produced the observed number of families with 
Z3pB = 200^-00 km. Instead, the required number of families 
with large parent bodies is systematically produced if the aster- 
oid belt was crossed by a large number of comets during the 
LHB, as expected in the Nice model (see Section |4|). However, 
for any reasonable size distribution of the cometary population, 
the same cometary flux that would produce the correct number of 
families with Dpb = 200^-00 km would produce too many fam- 
ilies with Dpb - 100 km relative to what is observed. Therefore, 
in the subsequent sections we look for mechanisms that might 
prevent detection of most of these families. 

More specifically, in Sec. [5] we discuss the possibility that 
families with Dpb - 100 km are so numerous that they can- 
not be identified because they overlap with each other. In Sec. [6] 
we investigate their possible dispersal below detectability due 
to the Yarkovsky effect and chaotic diffusion. In Sec. [7] we dis- 
cuss the role of the physical lifetime of comets. In Sec. [8] we 
analyse the dispersal of families due to the changes in the or- 
bits of the giant planets expected in the Nice model. In Sec. [9] 
we consider the subsequent collisional comminution of the fam- 
ilies. Of all investigated processes, the last one seems to be the 
most promising for reducing the number of visible families with 
Dps — 100 km while not affecting the detectability of old fami- 
lies with D PB = 200-400 km. 



Finally, in Section \W\ we analyse a curious portion of the 
main belt, located in a narrow semi-major axis zone bounded by 
the 5:2 and 7:3 resonances with Jupiter. This zone is severely 
deficient in small asteroids compared to the other zones of the 
main belt. For the reasons explained in the section, we think that 
this zone best preserves the initial asteroid belt population, and 
therefore we call it the "pristine zone". We checked the num- 
ber of families in the pristine zone, their sizes, and ages and we 
found that they are consistent with the number expected in our 
model invoking a cometary bombardment at the LHB time and a 
subsequent collisional comminution and dispersion of the family 
members. 

The conclusions follow in SectionfTTI 

2. A list of known families 

Although several lists of families exist in the literature (Zappala 
et al. 1995, Nesvorny et al. 2005, Parker et al. 2008, Nesvorny 
2010), we are going to identify the families once again. The rea- 
son is that we seek an upper limit for the number of old families 
that may be significantly dispersed and depleted, while the pre- 
vious works often focussed on well-defined families. Moreover, 
we need to calculate several physical parameters of the fami- 
lies (such as the parent-body size, slopes of the size-frequency 
distribution, a dynamical age estimate if not available in the liter- 
ature) which are crucial for further modelling. Last but not least, 
we use more precise synthetic proper elements from the AstDyS 
database (Knezevic & Milani 2003, version Aug 2010) instead 
of semi-analytic ones. 

We employed a hierarchical clustering method (HCM, 
Zappala et al. 1995) for the initial identification of families in 
the proper element space (a p , e p , sin/ p ), but then we had to per- 
form a lot of manual operations, because i) we had to select a 
reasonable cut-off velocity v C utoff, usually such that the number 
of members N{v cuto s) increases relatively slowly with increas- 
ing Vcutoff- ii) The resulting family should also have a "reason- 
able" shape in the space of proper elements, which should some- 
how correspond to the local dynamical featuresQiii) We checked 
taxonomic types (colour indices from the Sloan DSS MOC cat- 
alogue version 4, Parker et al. 2008), which should be consistent 
among family members. We can recognise interlopers or over- 
lapping families this way. iv) Finally, the size-frequency distri- 
bution should exhibit one or two well-defined slopes, otherwise 
the cluster is considered uncertain. 

Our results are summarised in online Tables [TJ-|3]and the po- 
sitions of families within the main belt are plotted in Figure [TJ 
Our list is "optimistic", so that even not very prominent families 
are included hereQ 

There are, however, several potential problems we are aware 

of: 

1 . There may be inconsistencies among different lists of fam- 
ilies. For example, sometimes a clump may be regarded as 
a single family or as two separate families. This may be the 
case of: Padua and Lydia, Rafita and Cameron. 

2. To identify families we used synthetic proper elements, 
which are more precise than the semi-analytic ones. 
Sometimes the families look more regular (e.g., Teutonia) 

1 For example, the Eos family has a complicated but still reasonable 
shape, since it is determined by several intersecting high-order mean- 
motion or secular resonances, see Vokrouhlicky et al. (2006). 

2 On the other hand, we do not include all of the small and less- 
certain clumps in a high-inclination region as listed by Novakovic et al. 
(201 1). We do not focus on small or high-/ families in this paper. 
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Fig. 1. Asteroids from the synthetic AstDyS catalogue plotted in the proper semimajor axis a p vs proper eccentricity e p (top panels) 
and fl p vs proper inclination sin 7 p planes (bottom panels). We show the identified asteroid families (left panels) with the positions 
of the largest members indicated by red symbols, and also remaining background objects (right panels). The labels correspond to 
designations of the asteroid families that we focus on in this paper. There are still some structures consisting of small objects in the 
background population, visible only in the inclinations (bottom right panel). These "halos" may arise for two reasons: i) a family has 
no sharp boundary and its transition to the background is smooth, or ii) there are bodies escaping from the families due to long-term 
dynamical evolution. Nevertheless, we checked that these halo objects do not significantly affect our estimates of parent-body sizes. 



or more tightly clustered (Beagle) when we use the syn- 
thetic elements. This very choice may, however, affect re- 
sults substantially! A clear example is the Teutonia family, 
which also contains the big asteroid (5) Astraea if the syn- 
thetic proper elements are used, but not if the semi-analytic 
proper elements are used. This is due to the large differences 
between the semi-analytic and synthetic proper elements of 
(5) Astraea. Consequently, the physical properties of the two 
families differ considerably. We believe that the family de- 
fined from the synthetic elements is more reliable. 
3. Durda et al. (2007) often claim a larger size for the parent 
body (e.g., Themis, Meliboea, Maria, Eos, Gefion), because 
they try to match the SFD of larger bodies and the results of 
SPH experiments. This way they also account for small bod- 
ies that existed at the time of the disruption, but which do not 
exist today since they were lost due to collisional grinding 
and the Yarkovsky effect. We prefer to use Dpurda instead of 
the value Dpb estimated from the currently observed SFD. 
The geometric method of Tanga et al. (1999), which uses 
the sum of the diameters of the first and third largest family 



members as a first guess of the parent-body size, is essen- 
tially similar to our approach. 

A complete list of all families' 
members is available at our web site 
http : / /sirrah . tro j a . mff . cuni . cz/~mira/mp/f ams/ 
including supporting figures. 

2.1. A definition of the production function 

To compare observed families to simulations, we define a "pro- 
duction function" as the cumulative number N(>D) of fami- 
lies with parent-body size Dpb larger than a given D. The ob- 
served production function is shown in Figure |2] and it is worth 
noting that it is very shallow. The number of families with 
Dpb - 100 km is comparable to the number of families in the 
D PB = 200^00 km range. 

It is important to note that the observed production function 
is likely to be affected by biases (the family sample may not be 
complete, especially below Dpb < 100 km) and also by long- 
term collisional/dynamical evolution which may prevent a de- 



3 



M. Broz et al.: Constraining the cometary flux throug] 



s 

I 10 





all observed families — ■ — 
catastrophic disruptions —* — 
^ + Qtarget + 5/3 Qproject 





\ ; 


— **** — **n m >H^ry*' " **■ ■ hi-. 







10 100 1000 

D[km] 

Fig. 2. A production function (i.e. the cumulative number N(>D) 
of families with parent-body size Dpb larger than D) for all 
observed families (black) and families corresponding to catas- 
trophic disruptions (red), i.e. with largest remnant/parent body 
mass ratio lower than 0.5. We also plot a theoretical slope ac- 
cording to Eq. ([1), assuming ^target = -3.2 and ^project = -1-2, 
which correspond to the slopes of the main belt population in the 
range D = 100-200 km and D = 15-60 km, respectively. 



tection of old comminutioned/dispersed families today (Marzari 
et al. 1999). 

From the theoretical point of view, the slope q of the produc- 
tion function N(>D) oc D q should correspond to the cumulative 
slopes of the size-frequency distributions of the target and pro- 
jectile populations. It is easy to show3 that the relation is 

5 

q — 2 + ^target + ^ ^project ■ (1) 

Of course, real populations may have complicated SFDs, with 
different slopes in different ranges. Nevertheless, any popula- 
tions that have a steep SFD (e.g. ^target = ^project = -2.5) would 
inevitably produce a steep production function (q = -4.7). 

In the following analysis, we drop cratering events and 
discuss catastrophic disruptions only, i.e. families which have 
largest remnant/parent body mass ratio less than 0.5. The rea- 
son is that the same criterion LR/PB < 0.5 is used in colli - 
sional models. Moreover, cratering events were not yet systemat- 
ically explored by SPH simulations due to insufficient resolution 
(Durda et al. 2007). 

2.2. Methods for family age determination 

If there is no previous estimate of the age of a family, we used 
one of the following three dynamical methods to determine it: 
i) a simple (a p , H) analysis as in Nesvorny et al. (2005); ii) a C- 
parameter distribution fitting as introduced by Vokrouhlicky et 
al. (2006); iii) a full N-body simulation described e.g. in Broz et 
al. (2011). 

In the first approach, we assume zero initial velocities, 
and the current extent of the family is explained by the size- 
dependent Yarkovsky semimajor axis drift. This way we can ob- 
tain only an upper limit for the dynamical age, of course. We 
show an example for the Eos family in Figure [3] The extent of 
the family in the proper semimajor axis vs the absolute magni- 
tude (fl p , H) plane can be described by the parametric relation 

0.2H = log 10 ^— - , (2) 



3 Assuming that the strength is approximately Q* D oc D 2 in the gravity 
regime, the necessary projectile size d oc (Q^) [I3 D (Bottke et al. 2005), 
and the number of disruptions n oc £) 2 £)i'mad' , P"'> ca . 



the asteroid belt during the late heavy bombardment 

Table 4. Nominal thermal parameters for S and C/X taxonomic 
types of asteroids: pbuik denotes the bulk density, p sur f the sur- 
face density, K the thermal conductivity, C± the specific thermal 
capacity, Aaond the Bond albedo and e the infrared emissivity. 



type 


Pbulk 


Psurf 


K 


C t h 


^Bond 


6 




(kg/m 3 ) 


(kg/m 3 ) 


(W/m/K) 


(J/kg/K) 






S 


2500 


1500 


0.001 


680 


0.1 


0.9 


C/X 


1300 


1300 


0.01 


680 


0.02 


0.9 



where a c denotes the centre of the family, and C is the parameter. 
Such relation can be naturally expected when the semimajor-axis 
drift rate is inversely proportional to the size, da/dt oc l/D, and 
the size is related to the absolute magnitude via the Pogson equa- 
tion H = -2.5 log 10 (/?i/D 2 /Dy), where Do denotes the reference 
diameter and py the geometric albedo (see Vokrouhlicky et al. 
2006 for a detailed discussion). The limiting value, for which 
all Eos family members (except interlopers) are above the cor- 
responding curve, is C = 1.5 to 2.0 x 10~ 4 AU. Assuming rea- 
sonable thermal parameters (summarised in Table 0, we calcu- 
late the expected Yarkovsky drift rates da/dt (using the theory 
from Broz 2006) and consequently can determine the age to be 
t < 1.5 to2.0Gyr. 

The second method uses a histogram N{C, C + AC) of the 
number of asteroids with respect to the C parameter defined 
above, which is fitted by a dynamical model of the initial ve- 
locity field and the Yarkovsky/YORP evolution. This enables us 
to determine the lower limit for the age too (so the resulting age 
estimate is t = 1 -S^ 5 Gyr for the Eos family). 

In the third case, we start an N-body simulation using a mod- 
ified SWIFT integrator (Levison and Duncan 1994), with the 
Yarkovsky/YORP acceleration included, and evolve a synthetic 
family up to 4 Gyr. We try to match the shape of the observed 
family in all three proper orbital elements (a p , e p , sin 7 p ). In prin- 
ciple, this method may provide a somewhat independent esti- 
mate of the age. For example, there is a 'halo' of asteroids in the 
surroundings of the nominal Eos family, which are of the same 
taxonomic type K, and we may fit the ratio Nha\o/N COK of the 
number of objects in the 'halo' and in the family 'core' (Broz et 
al., in preparation). 

The major source of uncertainty in all methods are unknown 
bulk densities of asteroids (although we use the most likely val- 
ues for the S or C/X taxonomic classes, Carry 2012). The age 
scales approximately as t oc p^ik- Nevertheless, we are still able 
to distinguish families that are young from those that are old, be- 
cause the allowed ranges of densities for S-types (2 to 3 g/cm 3 ) 
and C/X-types (1 to 2 g/cm 3 ) are limited (Carry 2012) and so are 
the allowed ages of families. 

2.3. Which families can be of LHB origin? 

The ages of the observed families and their parent-body sizes 
are shown in Figure |4] Because the ages are generally very un- 
certain, we consider that any family whose nominal age is older 
than 2 Gyr is potentially a family formed ~4Gyr ago, i.e. at the 
LHB time. If we compare the number of "young" (<2Gyr) and 
old families (>2Gyr) with Dpb = 200-400 km, we cannot see 
a significant over-abundance of old family formation events. On 
the other hand, we almost do not find any small old families. 

Only 12 families from the whole list may be possibly dated 
back to the late heavy bombardment, because their dynamical 
ages approach ~ 3.8 Gyr (including the relatively large uncer- 
tainties; see Table [5] which is an excerpt from Tables [TJ-[3]l. 
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AC = 0.1 x10" 4 AU 
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Fig. 3. An example of the Eos asteroid family, shown on the 
proper semimajor axis a p vs absolute magnitude H plot. We 
also plot curves defined by equation © and parameters a c = 
3.019 AU, C = 1.5 to 2.0 x 10~ 4 AU, which is related to the up- 
per limit of the dynamical age of the family. 
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Fig. 4. The relation between dynamical ages of families and the 
sizes of their parent bodies. Red labels correspond to catas- 
trophic disruptions, while cratering events are labelled in black. 
Some of the families are denoted by the designation of the 
largest member. The uncertainties of both parameters are listed 
in Tables QJ[3] (we do not include overlapping error bars here for 
clarity). 

If we drop cratering events and the families of Camilla 
and Hermione, which do not exist any more today (their exis- 
tence was inferred from the satellite systems, Vokrouhlicky et 
al. 2010), we end up with only five families created by catas- 
trophic disruptions that may potentially date from the LHB time 
(i.e. their nominal age is more than 2 Gy). As we shall see in 
Section|U this is an unexpectedly low number. 

Moreover, it is really intriguing that most "possibly-LHB" 
families are larger than Dpb - 200 km. It seems that old fam- 
ilies with Dpb - 100 km are missing in the observed sample. 
This is an important aspect that we have to explain, because it 
contradicts our expectation of a steep production function. 

3. Collisions in the main belt alone 

Before we proceed to scenarios involving the LHB, we try to 
explain the observed families with ages spanning 0^4 Gyr as a 
result of collisions only among main-belt bodies. To this pur- 
pose, we used the collisional code called Boulder (Morbidelli 
et al. 2009) with the following setup: the intrinsic probabilities 
Pi = 3.1 x 10~ 18 kirT 2 yr~', and the mutual velocities Vi mp = 
5.28 km/s for the MB vs MB collisions (both were taken from 
the work of Dahlgren 1998). The assumption of a single Vi mp 



Table 5. Old families with ages possibly approaching the LHB. 
They are sorted according to the parent body size, where Domda 
determined by the Durda et al. (2007) method is preferred to 
the estimate Dpb inferred from the observed SFD. An additional 
'c' letter indicates that we extrapolated the SFD down to D = 
km to account for small (unobserved) asteroids, an exclamation 
mark denotes a significant mismatch between Z?pb and £>Durda- 



des 


gnation 


D PB 




note 






(km) 


(km) 




24 


Themis 


209c 


380^130! 




10 


Hygiea 


410 


442 


cratering 


15 


Eunomia 


259 


292 


cratering 


702 


Alauda 


218c 


290-330! 


high-/ 


87 


Sylvia 


261 


272 


cratering 


137 


Meliboea 


174c 


240-290! 




375 


Ursula 


198 


240-280 


cratering 


107 


Camilla 


>226 




non-existent 


121 


Hermione 


>209 




non-existent 


158 


Koronis 


122c 


170-180 




709 


Fringilla 


99c 


130-140 


cratering 


170 


Maria 


100c 


120-130 





value is a simplification, but about 90 % collisions have mutual 
velocities between 2 and 8 km/s (Dahlgren 1998), which assures 
a similar collisional regime. 

The scaling law is described by the polynomial relation (r de- 
notes radius in cm) 

Gfl(r) = — (Qor a + Bpr") (3) 

with the parameters corresponding to basaltic material at 5 km/s 
(Benz & Asphaug 1999): 



p 


Go 


a B 


b 


<7fact 


(g/cm 3 ) 


(erg/g) 


(erg/g) 






3.0 


7x 10 7 


-0.45 2.1 


1.19 


1.0 



Even though not all asteroids are basaltic, we use the scaling law 
above as a mean one for the main-belt population. Below, we 
discuss also the case of significantly lower strengths (i.e. higher 
<7fact values). 

We selected the time span of the simulation 4 Gyr (not 
4.5 Gyr) since we are interested in this last evolutionary phase 
of the main belt, when its population and collisional activity is 
nearly same as today (Bottke et al. 2005). The outcome of a sin- 
gle simulation also depends on the "seed" value of the random- 
number generator that is used in the Boulder code to decide 
whether a collision with a fractional probability actually occurs 
or not in a given time step. We thus have to run multiple simula- 
tions (usually 100) to obtain information on this stochasticity of 
the collisional evolution process. 

The initial SFD of the main belt population conditions was 
approximated by a three-segment power law (see also thin grey 
line in Figure [5] 1st row) with differential slopes q a = -4.3 (for 
D > D\), qb = -2.2, q c - -3.5 (for D < D2) where the size 
ranges were delimited by D\ — 80 km and D2 = 16 km. We also 
added a few big bodies to reflect the observed shape of the SFD 
at large sizes (D > 400 km). The normalisation was N norm (D > 
D\) — 350 bodies in this case. 

We used the observed SFD of the main belt as the first con- 
straint for our collisional model. We verified that the outcome 
our model after 4 Gyr is not sensitive to the value of q c . Namely, 
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a change of q c by as much as + 1 does not affect the final SFD in 
any significant way. On the other hand, the values of the remain- 
ing parameters (q a , qt, D\, D2, N noim ) are enforced by the ob- 
served SFD. To obtain a reasonable fit, they cannot differ much 
(by more than 5-10 %) from the values presented above. 

We do not use only a single number to describe the num- 
ber of observed families (e.g. N - 20 for Dpb > 100 km), but 
we discuss a complete production function instead. The results 
in terms of the production function are shown in Figure |5] (left 
column, 2nd row). On average, the synthetic production func- 
tion is steeper and below the observed one, even though there is 
approximately a 5 % chance that a single realization of the com- 
puter model will resemble the observations quite well. This also 
holds for the distribution of Dpb = 200^400 km families in the 
course of time (age). 

In this case, the synthetic production function of Dpb > 
100 km families is not significantly affected by comminution. 
According to Bottke et al. (2005), most of D > 10 km fragments 
survive intact and a Dpb > 100 km family should be recognis- 
able today. This is also confirmed by calculations with Boulder 
(see Figure|5] left column, 3rd row). 

To improve the match between the synthetic and the ob- 
served production function, we can do the following: i) mod- 
ify the scaling law, or ii) account for a dynamical decay of the 
MB population. Using a substantially lower strength (qf dct = 5 
in Eq. (O, which is not likely, though) one can obtain a synthetic 
production function which is on average consistent with the ob- 
servations in the Dpb = 200^-00 km range. 

Regarding the dynamical decay, Minton & Malhotra (2010) 
suggest that initially the MB was three times more populous than 
today while the decay timescale was very short: after 100 Myr 
of evolution the number of bodies is almost at the current level. 
In this brief period of time, about 50 % more families will be 
created, but all of them will be old, of course. For the remaining 
~ 3.9 Gyr, the above model (without any dynamical decay) is 
valid. 

To conclude, it is possible - though not very likely - that 
the observed families were produced by the collisional activity 
in the main belt alone. A dynamical decay of the MB population 
would create more families that are old, but technically speaking, 
this cannot be distinguished from the LHB scenario, which is 
discussed next. 



4. Collisions between a "classical" cometary disk 
and the main belt 

In this section, we construct a collisional model and estimate an 
expected number of families created during the LHB due to col- 
lisions between cometary-disk bodies and main-belt asteroids. 
We start with a simple stationary model and we confirm the re- 
sults using a more sophisticated Boulder code (Morbidelli et al. 
2009). 

Using the data from Vokrouhlicky et al. (2008) for a "clas- 
sical" cometary disk, we can estimate the intrinsic collisional 
probability and the collisional velocity between comets and as- 
teroids. A typical time-dependent evolution of Pi and Vi mp is 
shown in Figure [6] The probabilities increase at first, as the 
transneptunian cometary disk starts to decay, reaching up to 
6 x 10~ 21 km -2 yr~', and after 100 Myr they decrease to zero. 
These results do not differ significantly from run to run. 




100 



f/ Myr 



Fig. 6. The temporal evolution of the intrinsic collisional proba- 
bility Pi (bottom) and mean collisional velocity Vi mp (top) com- 
puted for collisions between cometary-disk bodies and the main- 
belt asteroids. The time t = is arbitrary here; the sudden in- 
crease in Pi values corresponds to the beginning of the LHB. 

4.1. Simple stationary model 

In a stationary collisional model, we choose an SFD for the 
cometary disk, we assume a current population of the main belt; 
estimate the projectile size needed to disrupt a given target ac- 
cording to (Bottke et al. 2005) 



(2e;/v 2 mp ) 1/3 A a 



* disrupt 



(4) 



where Q* D denotes the specific energy for disruption and disper- 
sion of the target (Benz & Asphaug 1999); and finally calculate 
the number of events during the LHB as 



D 2 

^target 



^ target ^ Pi 



(t) «project(0 df ' 



(5) 



where n ta rget and n pro ject are the number of targets (i.e. main belt 
asteroids) and the number of projectiles (comets), respectively. 
The actual number of bodies (27,000) in the dynamical simula- 
tion of Vokrouhlicky et al. (2008) changes in the course of time, 
and it was scaled such that it was initially equal to the number of 
projectiles N(>ddi sm pt) inferred from the SFD of the disk. This 
is clearly a lower limit for the number of families created, since 
the main belt was definitely more populous in the past. 

The average impact velocity is VWj — lOkm/s, so we need 
the following projectile sizes to disrupt given target sizes: 



D target 


targets 




^disrupt for I s - = 3 to 6 

r /-'project 


(km) 


in the MB 


(J/kg) 


(km) 


100 


-192 


1 x 10 5 


12.6 to 23 


200 


-23 


4x 10 5 


40.0 to 73 



We tried to use various SFDs for the cometary disk (i.e., with 
various differential slopes q\ for D > Do and q2 for D < Do, 
the elbow diameter Do and total mass Mdi s k), including rather 
extreme cases (see Figure |7). The resulting numbers of LHB 
families are summarised in Table [6] Usually, we obtain sev- 
eral families with Dpb - 200 km and about 100 families with 
Dpb - 100 km. This result is robust with respect to the slope 
q2, because even very shallow SFDs should produce a lot of 
these familiesQ The only way to decrease the number of fam- 
ilies significantly is to assume the elbow at a larger diame- 
ter D - 150 km. 



4 The extreme case with q 2 = is not likely at all, e.g. because of 
the continuous SFD of basins on Iapetus and Rhea, which only ex- 
hibits a mild depletion of D a* 100 km size craters; see Kirchoff & 
Schenk (2010). On the other hand, Sheppard & Trujillo (2010) report 
an extremely shallow cumulative SFD of Neptune Trojans that is akin 
to low q2- 
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Fig. 5. Results of three different collisional models: main-belt alone which is discussed in Section [3] (left column), main-belt and 
comets from Section [4] (middle column), main-belt and disrupting comets from Section [7] (right column). 1st row: the initial and 
evolved size-frequency distributions of the main belt populations for 100 Boulder simulations; 2nd row: the resulting family pro- 
duction functions (in order to distinguish 100 lines we plot them using different colours ranging from black to yellow) and their 
comparison to the observations; 3rd row: the production function affected by comminution for a selected simulation; and 4th row: 
the distribution of synthetic families with Dpb > 50 km in the (age, Dpb) plot for a selected simulation, without comminution. The 
positions of synthetic families in the 4th-row figures may differ significantly for a different Boulder simulation due to stochasticity 
and low-number statistics. Moreover, in the middle and right columns, many families were created during the LHB, so there are 
many overlapping crosses close to 4 Gyr. 



It is thus no problem to explain the existence of approxi- 
mately five large families with Dpb = 200-400 km, which are 
indeed observed, since they can be readily produced during the 
LHB. On the other hand, the high number of Dpb - 100 km 
families clearly contradicts the observations, since we observe 
almost no LHB families of this size. 



4.2. Constraints from (4) Vesta 

The asteroid (4) Vesta presents a significant constraint for col- 
lisional models, being a differentiated body with a preserved 
basaltic crust (Keil 2002) and a 500 km large basin on its surface 
(a feature indicated by the photometric analysis of Cellino et al. 



1987), which is significantly younger than 4 Gyr (Marchi et al. 
2012). It is highly unlikely that Vesta experienced a catastrophic 
disruption in the past, and even large cratering events were lim- 
ited. We thus have to check the number of collisions between 
one D = 530km target and D ^ 35 km projectiles, which are 
capable of producing the basin and the Vesta family (Thomas et 
al. 1997). According to Table [6] the predicted number of such 
events does not exceed ~ 2, so given the stochasticity of the re- 
sults there is a significant chance that Vesta indeed experienced 
zero such impacts during the LHB. 
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Table 6. Results of a stationary collisional model between the cometary disk and the main belt. The parameters characterise the 
SFD of the disk: q\, q2 are differential slopes for the diameters larger/smaller than the elbow diameter Do, denotes the total 
mass of the disk, and « eve nts is the resulting number of families created during the LHB for a given parent body size Dps. The ranges 
of Events are due to variable density ratios ptai-get/Pproject = 1 to 3/1. 
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Fig. 7. Cumulative size-frequency distributions of the cometary 
disk tested in this work. All the parameters of our nominal choice 
are given in the top label; the other labels just report the param- 
eters that changed relative to our nominal choice. 

4.3. Simulations with the Boulder code 

To confirm results of the simple stationary model, we also per- 
formed simulations with the Boulder code. We modified the code 
to include a time-dependent collisional probability P;(f) and im- 
pact velocity V; mp (f) of the cometary-disk population. 

We started a simulation with a setup for the cometary disk 
resembling the nominal case from Table [6] The scaling law is 
described by Eq. (0, with the following parameters (the first set 
corresponds to basaltic material at 5 km/s, the second one to wa- 
ter ice, Benz & Asphaug 1999): 
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The intrinsic probabilities P\ = 3.1 x 10~ 18 krrT 2 yr~' and 
velocities Vi mp = 5.28 km/s for the MB vs MB collisions were 
again taken from the work of Dahlgren (1998). We do not ac- 
count for comet-comet collisions since their evolution is domi- 
nated by the dynamical decay. The initial SFD of the main belt 



was similar to the one in Section [3] q a = -4.2, qt, = -2.2, 
q c = -3.5, D[ = 80km, Di = 14km, and only the normalisation 
was increased up to N norm (D > D\) = 560 in this case. 

The resulting size-frequency distributions of 100 indepen- 
dent simulations with different random seeds are shown in 
Figure [5] (middle column). The number of LHB families (ap- 
proximately 10 with Dpb ^ 200 km and 200 with Dpb - 100 km) 
is even larger compared to the stationary model, as expected, be- 
cause we had to start with a larger main belt to get a good fit of 
the currently observed MB after 4 Gyr of collisional evolution. 

To conclude, the stationary model and the Boulder code give 
results that are compatible with each other, but that clearly con- 
tradict the observed production function of families. In partic- 
ular, they predict far too many families with D - 100 km par- 
ent bodies. At first sight, this may be interpreted as proof that 
there was no cometary LHB on the asteroids. Before jump- 
ing to this conclusion, however, one has to investigate whether 
there are biases against identifying of Dpb = 100 km families. 
In Sections |5]-[9] we discuss several mechanisms that all con- 
tribute, at some level, to reducing the number of observable 
Dpb = 100 km families over time. They are addressed in order 
of relevance, from the least to the most effective. 

5. Families overlap 

Because the number of expected Dpb > 100 km LHB families 
is very high (of the order of 100) we now want to verify if these 
families can overlap in such a way that they cannot be distin- 
guished from each other and from the background. We thus took 
192 main-belt bodies with D > 100 km and selected randomly 
100 of them that will break apart. For each one we created an ar- 
tificial family with 10 2 members, assume a size-dependent ejec- 
tion velocity V oc l/D (with V = 50m/s for D = 5 km) and 
the size distribution resembling that of the Koronis family. The 
choice of the true anomaly and the argument of perihelion at the 
instant of the break-up event was random. We then calculated 
proper elements (a p , e p , sin/ p ) for all bodies. This type of anal- 
ysis is in some respects similar to the work of Bendjoya et al. 
(1993). 

According to the resulting Figure [8] the answer to the ques- 
tion is simple: the families do not overlap sufficiently, and they 
cannot be hidden that way. Moreover, if we take only bigger bod- 
ies (D > 10 km), these would be clustered even more tightly. 
The same is true for proper inclinations, which are usually more 
clustered than eccentricities, so families could be more easily 
recognised. 
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Fig. 8. The proper semimajor axis a p vs the proper eccentric- 
ity e p for 100 synthetic families created in the main belt. It is the 
initial state, shortly after disruption events. We assume the size- 
frequency distribution of bodies in each synthetic family similar 
to that of the Koronis family (down to D ^ 2 km). Break-ups 
with the true anomaly / ^ to 30° and 150° to 180° are more 
easily visible on this plot, even though the choice of both / and 
the argument of perihelion m was random for all families. 

6. Dispersion of families by the Yarkovsky drift 

In this section, we model long-term evolution of synthetic fam- 
ilies driven by the Yarkovsky effect and chaotic diffusion. For 
one synthetic family located in the outer belt, we have performed 
a full N-body integration with the SWIFT package (Levison & 
Duncan 1994), which includes also an implementation of the 
Yarkovsky/YORP effect (Broz 2006) and second-order integra- 
tor by Laskar & Robutel (2001). We included 4 giant planets in 
this simulation. To speed-up the integration, we used ten times 
smaller sizes of the test particles and thus a ten times shorter 
time span (400 Myr instead of 4Gyr). The selected time step is 
Af = 91 d. We computed proper elements, namely their differ- 
ences Aflp, Ae p , A sin I p between the initial and final positions. 

Then we used a simple Monte-Carlo approach for the 
whole set of 100 synthetic families — we assigned a suitable 
drift Aflp(D) in semimajor axis, and also drifts in eccentric- 
ity A<?p and inclination A sin I p to each member of 100 families, 
respecting asteroid sizes, of course. This way we account for the 
Yarkovsky semimajor axis drift and also for interactions with 
mean-motion and secular resonances. This Monte-Carlo method 
tends to smear all structures, so we can regard our results as the 
upper limits for dispersion of families. 

While the eccentricities of small asteroids (down to D ^ 
2 km) seem to be dispersed enough to hide the families, there 
are still some persistent structures in inclinations, which would 
be observable today. Moreover, large asteroids (D > 10 km) 
seem to be clustered even after 4 Gyr, so that more than 50 % 
of families can be easily recognised against the background (see 
Figure|9j. We thus can conclude that it is not possible to disperse 
the families by the Yarkovsky effect alone. 

7. Reduced physical lifetime of comets in the MB 
crossing zone 

To illustrate the effects that the physical disruption of comets 
(due to volatile pressure build-up, amorphous/crystalline phase 
transitions, spin-up by jets, etc.) can have on the collisional 
evolution of the asteroid belt, we adopted here a simplistic as- 
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Fig. 9. The proper semimajor axis a p vs the proper inclination 
sin/p for 100 synthetic asteroid families (black dots), evolved 
over 4 Gyr using a Monte-Carlo model. The assumed SFDs cor- 
respond to the Koronis family, but we show only D > 10 km 
bodies here. We also include D > 10 km background asteroids 
(grey dots) for comparison. 



sumption. We considered that no comet disrupt beyond 1 .5 AU, 
whereas all comets disrupt the first time that they penetrate in- 
side 1.5 AU. Both conditions are clearly not true in reality: some 
comets are observed to blow up beyond 1.5 AU, and others are 
seen to survive on an Earth-crossing orbit. Thus we adopted our 
disruption law just as an example of a drastic reduction of the 
number of comets with small perihelion distance, as required to 
explain the absence of evidence for a cometary bombardment on 
the Moon. 

We then removed all those objects from output of comet evo- 
lution during the LHB that had a passage within 1.5 AU from 
the Sun, from the time of their first passage below this thresh- 
old. We then recomputed the mean intrinsic collision probabil- 
ity of a comet with the asteroid belt. The result is a factor ~3 
smaller than when no physical disruption of comets is taken into 
account as in Fig. [6] The mean impact velocity with asteroids 
also decreases, from 12km/s to 8 km/s. 

The resulting number of asteroid disruption events is thus de- 
creased by a factor ~4.5, which can be also seen in the produc- 
tion function shown in Figure [5] (right column). The production 
of families with Dpb = 200-400 km is consistent with observa- 
tions, while the number of Dpb - 100 km families is reduced to 
30-70, but is still too high, by a factor 2-3. More importantly, 
the slope of the production function remains steeper than that 
of the observed function. Thus, our conclusion is that physical 
disruptions of comets alone cannot explain the observation, but 
may be an important factor to keep in mind for reconciling the 
model with the data. 



8. Perturbation of families by migrating planets 
(a jumping-Jupiter scenario) 

In principle, families created during the LHB may be perturbed 
by still-migrating planets. It is an open question what the ex- 
act orbital evolution of planets was at that time. Nevertheless, a 
plausible scenario called a "jumping Jupiter" was presented by 
Morbidelli et al. (2010). It explains major features of the main 
belt (namely the paucity of high-inclination asteroids above the 
V6 secular resonance), and is consistent with amplitudes of the 
secular frequencies of both giant and terrestrial planets and also 
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Fig. 10. The proper semimajor axis vs the proper inclination for four synthetic families (distinguished by symbols) as perturbed by 
giant-planet migration. Left panel: the case when families were evolved over the "jump" due to the encounter between Jupiter and 
Neptune. Right panel: the families created just after the jump and perturbed only by later phases of migration. 



with other features of the solar system. In this work, we thus 
investigated this particular migration scenario. 

We used the data from Morbidelli et al. (2010) for the orbital 
evolution of giant planets. We then employed a modified SWIFT 
integrator, which read orbital elements for planets from an in- 
put file and calculated only the evolution of test particles. Four 
synthetic families located in the inner/middle/outer belt were in- 
tegrated. We started the evolution of planets at various times, 
ranging from fo to (to + 4 Myr) and stopped the integration at 
(fo + 4 My), in order to test the perturbation on families created 
in different phases of migration. Finally, we calculated proper 
elements of asteroids when the planets do not migrate anymore. 
(We also had to move planets smoothly to their exact current 
orbital positions.) 

The results are shown in Figure [10] While the proper eccen- 
tricities seem to be sufficiently perturbed and families are dis- 
persed even when created at late phases of migration, the proper 
inclinations are not very dispersed, except for families in the 
outer asteroid belt that formed at the very beginning of the gi- 
ant planet instability (which may be unlikely, as there must be a 
delay between the onset of planet instability and the beginning 
of the cometary flux through the asteroid belt). In most cases, the 
LHB families could still be identified as clumps in semi-major 
axis vs inclination space. We do not see any of such (a p , sin/ p )- 
clumps, dispersed in eccentricity, in the asteroid belt0 

The conclusion is clear: it is not possible to destroy low- 
e and low-/ families by perturbations arising from giant-planet 
migration, at least in the case of the "jumping-Jupiter" scenario^ 

9. Collisional comminution of asteroid families 

We have already mentioned that the comminution is not suffi- 
cient to destroy a Dps = 100 km family in the current environ- 
ment of the main belt (Bottke et al. 2005). 

However, the situation in case of the LHB scenario is differ- 
ent. Both the large population of comets and the several-times 
larger main belt, which has to withstand the cometary bombard- 
ment, contribute to the enhanced comminution of the LHB fam- 

3 High-inclination families would be dispersed much more owing to 
the Kozai mechanism, because eccentricities that are sufficiently per- 
turbed exhibit oscillations coupled with inclinations. 

6 The currently non-existent families around (107) Camilla and 
(121) Hermione — inferred from the existence of their satellites — 
cannot be destroyed in the jumping-Jupiter scenario, unless the fami- 
lies were actually pre-LHB and had experienced the jump. 



ilies. To estimate the amount of comminution, we performed 
the following calculations: i) for a selected collisional simula- 
tion, whose production function is close to the average one, we 
recorded the SFDs of all synthetic families created in the course 
of time; ii) for each synthetic family, we restarted the simula- 
tion from the time fo when the family was crated until 4 Gyr and 
saved the final SFD, i.e. after the comminution. The results are 
shown in Figure [TT] 

It is now important to discuss criteria, which enable us to 
decide if the comminutioned synthetic family would indeed be 
observable or not. We use the following set of conditions: Dpb > 
50 km, Dlf > 10 km (largest fragment is the first or the second 
largest body, where the SFD becomes steep), LR/PB < 0.5 (i.e. a 
catastrophic disruption). Furthermore, we define A^ mem bers as the 
number of the remaining family members larger than observa- 
tional limit Ai m ; t - 2 km and use a condition iV mem b er s > 10. The 
latter number depends on the position of the family within the 
main belt, though. In the favourable "almost-empty" zone (be- 
tween a p = 2.825 and 2.955 AU), N mem b as ^ 10 may be valid, 
but in a populated part of the MB one would need A^ m embers ^100 
to detect the family. The size distributions of synthetic families 
selected this way resemble the observed SFDs of the main-belt 
families. 

According to Figure [5] (3rd row), where we can see the 
production functions after comminution for increasing values 
of Mnembers, families with Dpb = 200-400 km remain more 
prominent than Dpb - 100 km families simply because they 
contain much more members with D > 10 km that survive in- 
tact. Our conclusion is thus that comminution may explain the 
paucity of the observed Dpb - 100 km families. 

10. "Pristine zone" between the 5:2 and 7:3 
resonances 

We now focus on the zone between the 5:2 and 7:3 mean-motion 
resonances, with a p = 2.825 to 2.955 AU, which is not as popu- 
lated as the surrounding regions of the main belt (see Figure [TJ. 
This is a unique situation, because both bounding resonances are 
strong enough to prevent any asteroids from outside to enter this 
zone owing the Yarkovsky semimajor axis drift. Any family for- 
mation event in the surroundings has only a minor influence on 
this narrow region. It thus can be called "pristine zone" because 
it may resemble the belt prior to creation of big asteroid families. 

We identified nine previously unknown small families that 
are visible on the (e p ,sin/ p ) plot (see Figure \T%. They are 
confirmed by the SDSS colours and WISE albedos, too. 
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Fig. 11. Left panel: the size-frequency distributions of the observed asteroid families. Middle panel: SFDs of 378 distinct synthetic 
families created during one of the collisional simulations of the MB and comets. Initially, all synthetic SFDs are very steep, in 
agreement with SPH simulations (Durda et al. 2007). We plot only the SFDs that fulfil the following criteria: Dpb > 50 km, 
Z?lf ^ 10 km, LR/PB < 0.5 (i.e. catastrophic disruptions). Right panel: the evolved SFDs after comminution. Only a minority of 
families are observable now, since the number of remaining members larger than the observational limit Di; m i t 2 km is often much 
smaller than 100. The SFD that we use for the simulation in Section[10]is denoted by red. 



a = 2 82 ■ 2 97 AU 



0.25 



sym > 



0.1 



0.25 



0.15 0.2 0.25 0.3 



a=2.82-2.97 AU 



0.15 0.2 0.25 0.3 




-0.4 -0.2 0.2 0.4 
a' [mag] 



Fig. 12. The "pristine zone" of the main belt (a p = 2.825 to 2.955 AU) displayed on the proper eccentricity e p vs proper inclina- 
tion sin I p plot. Left panel: the sizes of symbols correspond to the sizes of asteroids, the families are denoted by designations. Right 
panel: a subset of bodies for which SDSS data are available; the colours of symbols correspond to the SDSS colour indices a* and 
i - z (Parker et al. 2008). 



Nevertheless, there is only one big and old family in this zone 
(Dps > 100 km), i.e. Koronis. 

That at most one LHB family (Koronis) is observed in the 
"pristine zone" can give us a simple probabilistic estimate for 
the maximum number of disruptions during the LHB. We take 
the 192 existing main-belt bodies which have D > 100 km and 
select randomly 100 of them that will break apart. We repeat this 
selection 1000 times and always count the number of families in 
the pristine zone. The resulting histogram is shown in Figure[T3l 
As we can see, there is very low (<0.001) probability that the 
number of families in the pristine zone is zero or one. On average 
we get eight families there, i.e. about half of the 16 asteroids 
with D > 100 km present in this zone. It seems that either the 
number of disruptions should be substantially lower than 100 or 
we expect to find at least some "remnants" of the LHB families 
here. 

It is interesting that the SFD of an old comminutioned family 
is very flat in the range D — 1 to 10 km (see Figure fTTb — simi- 
lar to those of some of the "less certain" observed families! We 



may speculate that the families like (918) Itha, (5567) Durisen, 
(12573) 1999 NJ 53 , or (15454) 1998 YB 3 (all from the pristine 
zone) are actually remnants of larger and older families, even 
though they are denoted as young. It may be that the age esti- 
mate based on the (a p , H) analysis is incorrect since small bodies 
were destroyed by comminution and spread by the Yarkovsky ef- 
fect too far away from the largest remnant, so they can no longer 
be identified with the family. 

Finally, we have to ask an important question: what does 
an old/comminutioned family with Dpb - 100 km look like in 
proper-element space? To this aim, we created a synthetic fam- 
ily in the "pristine zone", and assumed the family has A^ mem beis - 
100 larger than Airmt - 2 km and that the SFD is already flat in 
the D — 1 to 10 km range. We evolved the asteroids up to 4Gyr 
due to the Yarkovsky effect and gravitational perturbations, us- 
ing the N-body integrator as in Section [6] Most of the D 2 km 
bodies were lost in the course of the dynamical evolution, of 
course. The resulting family is shown in Figure [14] We can also 
imagine that this family is placed in the pristine zone among 
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Fig. 13. The histogram for the expected number of LHB families 
located in the "pristine zone" of the main belt. 
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Fig. 14. The proper eccentricity vs proper inclination of one 
synthetic old/comminutioned family evolved dynamically over 
4Gyr. Only a few family members (N - 10 1 ) remained from 
the original number of N(D > 2 km) =s 10 2 . The scales are the 
same as in Figure [12] so we can compare it easily to the "pris- 
tine zone". 

other observed families, to get a feeling of whether it is easily 
observed or not (refer to Figure [12}. 

It is clear that such family is hardly observable even in the 
almost empty zone of the main belt! Our conclusion is that the 
comminution (as given by the Boulder code) can explain the 
paucity of Dpb - 100 km LHB families, since we can hardly 
distinguish old families from the background. 

11. Conclusions 

In this paper we investigated the cometary bombardment of the 
asteroid belt at the time of the LHB, in the framework of the Nice 
model. There is much evidence of a high cometary flux through 
the giant planet region, but no strong evidence of a cometary 
bombardment on the Moon. This suggests that many comets 
broke up on their way to the inner solar system. By investigat- 
ing the collisional evolution of the asteroid belt and comparing 
the results to the collection of actual collisional families, our aim 
was to constrain whether the asteroid belt experienced an intense 
cometary bombardment at the time of the LHB and, if possible, 
constrain the intensity of this bombardment. 

Observations suggest that the number of collisional families 
is a very shallow function of parent-body size (that we call in 



this paper the "production function"). We show that the colli- 
sional activity of the asteroid belt as a closed system, i.e. without 
any external cometary bombardment, in general does not pro- 
duce such a shallow production function. Moreover, the number 
of families with parent bodies larger than 200 km in diameter 
is in general too small compared to the observations. However, 
there is a lot of stochasticity in the collisional evolution of the 
asteroid belt, and about 5 % of our simulations actually fit the 
observational constraints (shallowness of the production func- 
tion and number of large families) quite well. Thus, in principle, 
there is no need for a bombardment due to external agents (i.e. 
the comets) to explain the asteroid family collection, provided 
that the real collisional evolution of the main belt was a "lucky" 
one and not the "average" one. 

If one accounts for the bombardment provided by the comets 
crossing the main belt at the LHB time, predicted by the Nice 
model, one can easily justify the number of observed families 
with parent bodies larger than 200 km. However, the resulting 
production function is steep, and the number of families pro- 
duced by parent bodies of 100 km is almost an order of magni- 
tude too large. 

We have investigated several processes that may decimate 
the number of families identifiable today with 100 km parent 
bodies, without considerably affecting the survival of families 
formed from larger parent bodies. Of all these processes, the col- 
lisional comminution of the families and their dispersal by the 
Yarkovsky effect are the most effective ones. Provided that the 
physical disruption of comets due to activity reduced the effec- 
tive cometary flux through the belt by a factor of * 5, the result- 
ing distribution of families (and consequently the Nice model) is 
consistent with observations. 

To better quantify the effects of various cometary-disruption 
laws, we computed the numbers of asteroid families for dif- 
ferent critical perihelion distances q CI \t and for different disrup- 
tion probabilities p a -\ t of comets during a given time step (At = 
500 yr in our case). The results are summarised in Figure [15] 
Provided that comets are disrupted frequently enough, namely 
the critical perihelion distance has to be at least q cr it Z 1 AU, 
while the probability of disruption is p cl i t = 1, the number of 
Z)pB > 100 km families drops by the aforementioned factor of 
ss 5. Alternatively, q cr \ t may be larger, but then comets have to 
survive multiple perihelion passages (i.e. p cr \ t have to be lower 
than 1). It would be very useful to test these conditions by in- 
dependent models of the evolution and physical disruptions of 
comets. Such additional constraints on cometary-disruption laws 
would then enable study of the original size-frequency distribu- 
tion of the cometary disk in more detail. 

We can also think of two "alternative" explanations: i) phys- 
ical lifetime of comets was strongly size-dependent so that 
smaller bodies break up easily compared to bigger ones; ii) high- 
velocity collisions between hard targets (asteroids) and very 
weak projectiles (comets) may result in different outcomes than 
in low-velocity regimes explored so far. Our work thus may also 
serve as a motivation for further SPH simulations. 

We finally emphasize that any collisional/dynamical models 
of the main asteroid belt would benefit from the following ad- 
vances: 

i) determination of reliable masses of asteroids of various 
classes. This may be at least partly achieved by the Gaia 
mission in the near future. Using up-to-date sizes and shape 
models (volumes) of asteroids one can then derive their den- 
sities, which are directly related to ages of asteroid families. 
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Fig. 15. The numbers of collisional families for different critical 
perihelion distances q CI it at which comets break up and disrup- 
tion probabilities p CI \t during one time step (At = 500 yr). In the 
top panel, we vary q a i t while keeping p cr j t = 1 constant. In the 
bottom panel, q cr \ t = 1.5 AU is constant and we vary p a i V We al- 
ways show the number of catastrophic disruptions with parent- 
body sizes Dpb > 100 km (red line) and 200km (black line). 
The error bars indicate typical (1-cr) spreads of Boulder simu- 
lations with different random seeds. The observed numbers of 
corresponding families are indicated by thin dotted lines. 

ii) Development of methods for identifying asteroid families 
and possibly targeted observations of larger asteroids ad- 
dressing their membership, which is sometimes critical for 
constructing size-frequency distributions and for estimating 
parent-body sizes. 

iii) An extension of the SHP simulations for both smaller and 
larger targets, to assure that the scaling we use now is 
valid. Studies and laboratory measurements of equations of 
states for different materials (e.g. cometary-like, porous) are 
closely related to this issue. 

The topics outlined above seem to be the most urgent develop- 
ments to be pursued in the future. 
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Table 1. A list of asteroid families and their physical parameters. There are the following columns: v cu toff is the selected cut-off velocity for the hierarchical clustering, N the 
corresponding number of family members, pv the adopted value of the geometric albedo for family members which do not have measured diameters (from Tedesco et al. 2002 
or Masiero et al. 201 1, a letter 'w' indicates it was necessary to use the WISE data to obtain median/mean albedo), taxonomic classification (according to the Sloan DSS MOC 4 
colours, Parker et al. 2008), Dpb parent body size, an additional 'c' letter indicates that we prolonged the SFD slope down to zero D (a typical uncertainty is 10 %), Dom-da PB 
size inferred from SPH simulations (Durda et al. 2007), an exclamation mark denotes a significant mismatch with Dpb, LR/PB the ratio of the volumes of the largest remnant to 
the parent body (an uncertainty corresponds to the last figure, a range is given if both D PB and D Durda are known), v esc the escape velocity, q x the slope of the SFD for larger D, 
q2 the slope for smaller D (a typical uncertainty of the slopes is 0.2, if not indicated otherwise), dynamical age including its uncertainty. 
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174c 


24U-29U ! 


n cn n on 

0.59-0.20 


102 


-1.9 


-1.8 


<3.0 


old.' 


142 


Polana (Nysa) 


An 
Ol) 


1 A A Q 

3443 


n ncc, 
U.U55W 


L. 


/5 


if 


n /i o 
U.42 


A C 

45 


h n 
- /.U 


-3.6 


<1.5 


overlaps with Nysa 


145 


Adeona 


cn 

50 


1 1 A 1 

1 lot 


n nAC 

0.065 


c 


ni„ 

171c 


1 o c 


n An n c /i 

0.69-0.54 


1 m 
101 


c o 

-5.2 


-2.8 


0.7 + 0.5 


cut by J5/2 resonance 


158 


Koronis 


cn 

50 


4225 


n 1 /i "7 

0.147 


c 


1 oo~ 

122c 


17U-180 


n no /i n nnn 

0.024-0. 009 


AO 

68 


-3.6(0.3) 


-2.3 


2.5 ± 1.0 


LHB? 


163 


Erigone 


An 

60 


1059 


n ncA 

0.056 


c/x 


79 


114 


n 7n n ^ic. 
U. /9-(J.2D 


/1 A 

46 




-3.6 


0.3 ± 0.2 




170 


Maria 


on 


Qnn/i 
3U94 


n o /i n.,. 
U.249W 


c 


1U/C 


1 on 1 on 


n mn n n<i o 
U.U/U-U.U45 


AQ 

03 


o c/n i\ 
-2.5(11.3) 


-2.8 


3.0+1.0 


T UTJ ^ 


221 


Eos 


5U 


cm a 
59/0 


n i Qn 
U. IjU 


K 


onO/i 
2U8C 


ool ! 


n 1 q n mn 


1 OQ 

123 


1 c 

-3.5 


-2.1 


1.3 ± 0.2 




283 


Emma 


/5 


1 A C 

345 


n ncn 
U.U5U 




1 CO 

152 


1 


n no a 
U.y2-(J.51 


nn 
9U 




-3.2 


<1.0 


satellite 


293 


Brasilia 


An 
DU 


TOO 

282 


n 1 t?,„ 
U. 1 /5W 


C/A 


"2 A 
34 


^ 


n non 
U.U2U 


on 
2U 


i /i /"n c\ 
— 1.4(U.5) 


-3.7 


0.05 ± 0.04 


(293) is interloper 


363 


Padua (Lydia) 


51) 


cnA 
590 


n nm 
U.U9/ 


L./A 


"7 A 
10 


lUo 


n n/i c n n 1 t 
U.U4D-U.U1 / 


/I c 

45 


1 o 
— 1.0 


-3.0 


0.3 ± 0.2 




396 


Aeolia 


on 
2U 


124 


n i "7 1 
U.l 11 


(_/A 


1 C 

35 




n n^:<c n nr\ 
U.VOD-U. /U 


on 
2U 




-4.3 


<0.1 


cratering 


410 


Chloris 


nn 
9U 


ocn 

259 


n nc7 
U.U5 / 


/ - 
L 


1 OA„ 

120C 


1 ZA 


n nco n co 


T/1 

/4 




-2.1 


0.7 ± 0.4 


490 


Veritas 












1UU-1 / / 










0.0083 ± 0.0005 


//inn\ if i;i^/ii./ in^j/m^r Ayi:^ii/.i ^» „i om i \ 

(49U) is liKely interloper (Micnel et al. 2U1 1 ) 


569 


Misa 


;U 


C/1Q 

543 


n nQ i 
U.UJl 


L. 


8SC 


11/ 


n co n oc 


CO 

52 


1 n 

— 3.y 


-2.3 


0.5 + 0.2 




606 


Brangane 


30 


81 


0.102 


s 


37 


46 


n no r\ ao 


22 




-3.8 


0.05 ± 0.04 




668 


Dora 


cn 

50 


837 


0.054 


c 


85 


165! 


n m 1 n nn /i 

U.U31-U.UU4 


50 


-4.2 


-1.9 


0.5 ± 0.2 




808 


Merxia 


50 


549 


0.227 




37 


121 ! 


66-0 01 8 


22 


-2.7 


-3.4 


0.3 ± 0.2 




832 


Karin 








s 




40 










0.0058 ± 0.0002 




845 


Naema 


30 


173 


0.081 


c 


77c 


81 


0.35-0.30 


46 


-5.2 


-2.9 


0.1+0.05 




847 


Agnia 


40 


1077 


0.177 


s 


39 


61 


0.38-0.10 


23 


-2.8 


-3.1 


0.2 + 0.1 




1128 


Astrid 


50 


265 


0.079 


c 


43c 


? 


0.52 


25 


-1.7 


-2.6 


0.1+0.05 




1272 


Gefion 


60 


19477 


0.20 


s 


74c 


100-150! 


0.001-0.004 


60 


-4.3 


-2.5 


0.48 ± 0.05 


Nesvorny et al. (2009), L chondrites 


1400 


Tirela 


80 


1001 


0.070 


s 


86 




0.12 


86 


-4.2 


-3.4 


<1.0 




1658 


Innes 


70 


621 


0.246w 


s 


27 


? 


0.14 


16 


-4.9 


-3.5 


<0.7 


(1644) Rafita is interloper 


1726 


Horfrneister 


40 


822 


0.035 


c 


93c 


134 


0.022-0.007 


55 


-4.5 


-2.7 


0.3 ± 0.2 




3556 


Lixiaohua 


60 


439 


0.044w 


c/x 


62 


220! 


0.029-0.001 


35 


-6.1 


-3.3 


0.15 + 0.05 


Novakovic et al. (2010) 


3815 


Konig 


60 


177 


0.044 


c 


33 


? 


0.32 


20 


? 


-3.0 


<0.1 


(1639) Bower is interloper 


4652 


Iannini 








s 














0.005 ± 0.005 




9506 


Telramund 


40 


146 


0.217w 


s 


22 




0.05 


13 


-3.9 


-3.7 


<0.5 




18405 


1993 FY 12 


50 


44 


0.171w 


c/x 


15 




0.23 


15 


-2.4 


-2.4 


<0.2 


cut by J5/2 resonance 



Table 2. Continuation of Tabled] 



designation 


VcutofF 

m/s 


N 


Pv 


tax. 


km 


£>Dunia 

km 


LR/PB 


"esc 

rn/s 


?i 


?2 


age 
Gyr 


notes, references 


158 


Koronis ( 2) 








c 

"J 


j J 












n m s -t-0 nns 

U.U1J ± U.UUJ 


T~- on 5— jj 

cratcnng, JVlolnar & rlacgsrt (ZUU7) 


298 


Baptistina 


SO 


1 740 


n 1 fifw 

U. 1 WW 








n 17 

u. 1 / 


9 1 
z 1 


—3 fi 

J.U 


—7 4 


<-ft 3 

^U. J 


jjdiL ui Liic rujid iciiiiiiy 


434 


Hungaria 


7 on 
zuu 


4SQ8 


n 3S 

U.J J 




9S 




n 1 s 

U.1J 


1 s 


S Q 

J.y 


3 1 
J. 1 


n s -1- n 7 
U.J ± u.z 


WdlllCl CL dl . ^ZUIU^J 


627 


Charis 


80 

oU 


73S 

ZjJ 


n ns 1 

U.Uol 


c 







n S3 

U.J J 


^S 

jj 


7 


3 A 
J.H 


^ 1 n 

V 1 .u 




778 


1 11CU UdlUd 


8S 

OJ 


1 S4 
1 JH 


n nfin 

U.UOU 




7 /C 




n 9Q 


S7 


'7 


7 Q 


n nn7 -1- n nn7 
u.uu / ± u.uuz 


cratcnng, JNovskovic (zUlu) 


302 


Clarissa 


30 


75 


0.054 


c 


39 


- 


0.96 


23 


7 


-3.1 


<0.1 


cratering, Nesvomy (2010) 


656 


Beagle 


24 


63 


0.089 


c 


64 


- 


0.56 


38 


-1.3 


-1.4 


<0.2 




752 


Sulamitis 


fin 


1 1 


047 

U.UH-Z 


c 


fiS 




n 83 

U.o J 


J J 


— fi s 

U.J 


—7 3 

— Z.J 


^n 4 




1189 


Terentia 


50 


18 


0.070 




56 




0.990 


33 


7 


-2.6? 


<0.2 




1892 


Lucienne 


i on 


S7 


773w 


c 


1 t 




n 71 

u. / 1 






7 


—4 4 


^n 3 

^U. J 




7353 


Kazvia 


sn 

ju 


~>3 
_ J 


U.ZUDW 


c 



i fi 
1 




n S7 

U.J / 


Q 



9 


1 8 
1.0 


^n 1 
<u. 1 




10811 


Lau 


100 


15 


0.273w 


s 


11 




0.77 


5 


? 


-2.8 


<0.1 




18466 


1995 SU 37 


40 


71 


0.241w 


s 


14 


- 


0.045 


7 


? 


-5.0 


<0.3 




1270 


Datura 








c 














n nno4s n nnnfin 

U.UUUtJ u.uuuuu 


1UCI1L111CU 111 UaLUlaLlllg-ClClllCllL apaL-C, 


14627 


Emilkowalski 








CIY 














n nnm q n nnn7s 

u.uuu 1 y u.uuuz j 


INCftVUlliy OC VUKJULlllllLKy ^zuuu^ 


16598 


1992 YC 2 






















00005-0 00075 




21509 


Lucascavin 






















000^-0 0008 




2384 


Schulhof 








c 














0007 000Q 
u.uuu / -u.uuu? 


VURJUUllllClvy GL fHCSVUlliy ^ZUll^ 


27 


Euterpe 


70 


268 


0.260w 


s 


118c 


- 


0.998 


70 


-2.9 


-2.2 


<1.0 


cratering, Parker et al. (2008) 


375 


Ursula 


80 


111 


0.057w 


c 


203c 


240-280 


0.71-0.43 


120 


-4.1 


-2.3 


<3.5 


old? 


1044 


Teutonia 


so 




U. J'+J 


c 


97 1 90 
z / - 1ZU 




o 1 7 o qs 

U. 1 /-U.70 


1 f\ 1 1 

ID- / 1 


3 s 

j.j 


^ Q 


^-0 s 

<U. J 


depends on (5) Astraea membership 


1296 


Andree 


60 


H-U 1 


9QOw 


c 


17 74. 




010 QS 

U.U 1U-U.7 J 


10 zn 

1 U HO 


7 


—9 QfO 


<r"l 
1 .u 


ucpciius Ull \ /7) TjUiyilUlllC iiiciiiucisiiip 


2007 


McCuskey 


"XA 


9^6 


u.uo 


r 1 


90 




n zti 


1 7 

1 / 


7 


— J.D 


<-U. J 


overlaps with Nysa/Polana 


2085 


Henan 


J4 




u.zuu w 


c 


11 
z / 




n 1 1 


1 u 


—4. 7 


— J.Z 


<rl 

v. 1 .U 




2262 


Mitidika 


OJ 


H- 1U 


U.UOH-W 




H-7- 17^ 




n n^7 n s 1 


9A Af\ 


4. S 
H.J 


9 9 
— z.z 


< l.U 


depends on \ /oj) ^wetana membersnip, 




























^ZZUZ^ la UlLCllUpCl, UVCllaps W1L11 J LlllU 


2 


Pallas 


200 


64 


0.163 


B 


498c 


- 


0.9996 


295 


7 


-2.2 


<0.5 


high-/, Carruba (2010) 


25 


Phocaea 


160 


1370 


0.22 


S 


92 




0.54 


55 


-3.1 


-2.4 


<2.2 


old? high-//e, cut by v 6 resonance, Carruba (2009) 


148 


Gallia 


150 


S7 

j 1 


n 1 fiQ 

u. 1 \jy 


Q 


08 




n oss 

U.UJO 


S8 

JO 


7 


-3 fi 

J.U 


^n 4s 


high-/ 


480 


Hansa 


i sn 


fiS 1 
j 1 


f) 7Sfi 


e 


fin 

uu 




0.83 


j j 


—4 Q 


— 3 7 

— J.Z 


^•1 fi 
V. 1 .u 


hi oh I 

ingii i 


686 


Gersuind 


130 


178 


0.146 


S 


52c 




0.48 


40 


7 


-2.7 


<0.8 


high-/, Gil-Hutton (2006) 


945 


Barcelona 


110 


129 


0.248 


s 


28 




0.77 


16 


? 


-3.5 


<0.35 


high-/, Foglia & Masi (2004) 


1222 


Tina 


110 


37 


0.338 


s 


21 




0.94 


12 


7 


-4.1 


<0.15 


high-/ 


4203 


Brucato 






















<1.3 


in freq. space 


31 


Euphrosyne 


100 


851 


0.056 


c 


259 




0.97 


153 


-4.9 


-3.9 


<1.5 


cratering, high-/, Foglia & Massi (2004) 


702 


Alauda 


120 


791 


0.070 


B 


218c 


290-330! 0.025 


129 


-3.9 


-2.4 


<3.5 


old? high-/, cut by J2/1 resonance, satellite 




























(Margot & Rojo 2007) 


107 


Camilla 


? 


? 


0.054 




>226 


? 


? 


7 


? 


7 


3.8? 


LHB? Cybele region, non-existent today, 


121 


Hermione 


9 


? 


0.058 




>209 


7 


? 


7 


7 


7 


3.8? 


LHB ? Vokrouhlicky et al. (20 1 0) 



Table 3. Continuation of Table Q] 



desi 


gnation 


^cutoff 


N 


Pv 


tax. 


D PB 


^Durda 


LR/PB 


Vesc 


?i 


12 


age 


notes, references 






m/s 








km 


km 




m/s 






Gyr 




1303 


Luthera 


100 


142 


0.043 


X 


92 




0.81 


54 


-3.9 


-2.7 


<0.5 


above (375) Ursula 


1547 


Nele 


20 


57 


0.311w 


X 


19 




0.85 


11 


? 


-2.8(0.3) 


<0.04 


close to (3) Juno 


2732 


Witt 


60 


985 


0.260w 


s 


25 




0.082 


15 


-4.0(0.3) 


-3.8 


<1.0 


only part with sin / > 0.099, above (363) Padua 


81 


Terpsichore 


120 


70 


0.052 


c 


119 




0.993 


71 


? 


-4.4 


<0.5 


cratering, less-certain families in the "pristine zone" 


709 


Fringilla 


140 


60 


0.047 


X 


99c 


130-140 


0.93-0.41 


59 


-6.2 


-1.7 


<2.5 


old? 


918 


Itha 


140 


63 


0.23 


s 


38 




0.16 


22 


-2.7 


-1.5 


<1.5 


shallow SFD 


5567 


Durisen 


100 


18 


0.044w 


X 


42 




0.89 


25 


7 


-1.7 


<0.5 


shallow SFD 


5614 


Yakovlev 


100 


34 


0.05 


c 


22 




0.28 


13 


? 


-3.2 


<0.2 




12573 


1999 Nfc 


40 


13 


0.1 90w 


c 


15 




0.13 


9 


? 


-2.0(0.5) 


<0.6 


incomplete SFD 


15454 


1998 YB 3 


50 


14 


0.054w 


c 


21 




0.41 


13 


? 


-1.6(0.3) 


<0.5 


shallow SFD 


15477 


1999 CGi 


110 


144 


0.098w 


s 


25 




0.065 


14 


? 


-4.6(0.5) 


<1.5 




36256 


1999 XT 17 


60 


30 


0.210w 


s 


17 




0.037 


10 


? 


-1.4(0.5) 


<0.3 


shallow SFD 



